Go back to the Convergence Rates page.
Let us set some global options for all code chunks in this document.
knitr::opts_chunk$set(
message = FALSE, # Disable messages printed by R code chunks
warning = FALSE, # Disable warnings printed by R code chunks
echo = TRUE, # Show R code within code chunks in output
include = TRUE, # Include both R code and its results in output
eval = TRUE, # Evaluate R code chunks
cache = FALSE, # Enable caching of R code chunks for faster rendering
fig.align = "center",
#out.width = "100%",
retina = 2,
error = TRUE,
collapse = FALSE
)
rm(list = ls())
set.seed(1982)# Install R-INLA package
# install.packages("INLA",repos = c(getOption("repos"),INLA ="https://inla.r-inla-download.org/R/testing"), dep = TRUE)
# Update R-INLA package
# inla.upgrade(testing = TRUE)
# Install inlabru package
# remotes::install_github("inlabru-org/inlabru", ref = "devel")
# Install rSPDE package
# remotes::install_github("davidbolin/rspde", ref = "devel")
# Install MetricGraph package
# remotes::install_github("davidbolin/metricgraph", ref = "devel")
library(INLA)
library(inlabru)
library(rSPDE)
library(MetricGraph)
library(Matrix)
library(dplyr)
library(plotly)
library(scales)
library(patchwork)
library(tidyr)
library(ggplot2)
library(reshape2)
library(sf)
library(here)
library(rmarkdown)
library(knitr)
library(grateful) # Cite all loaded packages
library(latex2exp)
library(plotrix)
rm(list = ls()) # Clear the workspace
set.seed(1982) # Set seed for reproducibility# Eigenfunctions for the tadpole graph
tadpole.eig <- function(k,graph){
x1 <- c(0,graph$get_edge_lengths()[1]*graph$mesh$PtE[graph$mesh$PtE[,1]==1,2])
x2 <- c(0,graph$get_edge_lengths()[2]*graph$mesh$PtE[graph$mesh$PtE[,1]==2,2])
if(k==0){
f.e1 <- rep(1,length(x1))
f.e2 <- rep(1,length(x2))
f1 = c(f.e1[1],f.e2[1],f.e1[-1], f.e2[-1])
f = list(phi=f1/sqrt(3))
} else {
f.e1 <- -2*sin(pi*k*1/2)*cos(pi*k*x1/2)
f.e2 <- sin(pi*k*x2/2)
f1 = c(f.e1[1],f.e2[1],f.e1[-1], f.e2[-1])
if((k %% 2)==1){
f = list(phi=f1/sqrt(3))
} else {
f.e1 <- (-1)^{k/2}*cos(pi*k*x1/2)
f.e2 <- cos(pi*k*x2/2)
f2 = c(f.e1[1],f.e2[1],f.e1[-1],f.e2[-1])
f <- list(phi=f1,psi=f2/sqrt(3/2))
}
}
return(f)
}
# Function to compute the true covariance matrix
gets_true_cov_mat <- function(graph, kappa, tau, alpha, n.overkill){
Sigma.kl <- matrix(0,nrow = dim(graph$mesh$V)[1],ncol = dim(graph$mesh$V)[1])
for(i in 0:n.overkill){
phi <- tadpole.eig(i,graph)$phi
Sigma.kl <- Sigma.kl + (1/(kappa^2 + (i*pi/2)^2)^(alpha))*phi%*%t(phi)
if(i>0 && (i %% 2)==0){
psi <- tadpole.eig(i,graph)$psi
Sigma.kl <- Sigma.kl + (1/(kappa^2 + (i*pi/2)^2)^(alpha))*psi%*%t(psi)
}
}
Sigma.kl <- Sigma.kl/tau^2
return(Sigma.kl)
}edge1 <- rbind(c(0,0), c(1,0))
theta <- seq(from = -pi, to = pi,length.out = 100)
edge2 <- cbind(1+1/pi+cos(theta)/pi, sin(theta)/pi)
edges <- list(edge1, edge2)
graph <- metric_graph$new(edges = edges)
graph$plot() +
ggtitle("Tadpole graph") +
theme_minimal() +
theme(text = element_text(family = "Palatino")) +
coord_fixed(ratio = 1)# parameters
h.ok <- 2^-9
type <- "covariance"
type_rational_approximation = "brasil"
rho <- 0.5
#m = 4
sigma <- 1
n.overkill = 1000
# Mesh sizes
h_aux <- c(6:3)
h_vector <- 2^-h_aux
h_label <- paste0("2^-", h_aux, "")
h_label_latex <- sprintf("$2^{-%d}$", h_aux)
# Beta values
beta_aux <- c(4, 4.5, 5, 5.5, 6, 6.5, 7)
beta_vector <- beta_aux/8
beta_label <- paste0(beta_aux, "/8")
theoretical_rate <- pmin(4*beta_vector-1/2,2)graph.ok <- graph$clone()
# Build graph with overkill mesh
graph.ok$build_mesh(h = h.ok)
graph.ok$plot(mesh = TRUE) +
ggtitle("Tadpole graph with overkill mesh") +
theme_minimal() +
theme(text = element_text(family = "Palatino")) +
coord_fixed(ratio = 1)Each A[[i]] below is the projection matrix corresponding
to the graph with mesh size h_vector[i] onto the overkill
mesh.
# Initialize the list of graphs and the list of projection matrices
graphs <- list()
A <- list()
for(i in 1:length(h_vector)){
graphs[[i]] <- graph$clone()
graphs[[i]]$build_mesh(h = h_vector[i])
A[[i]] <- graphs[[i]]$fem_basis(loc.ok)
}
# Print the dimensions of the projection matrices
print(lapply(A, dim))## [[1]]
## [1] 1536 192
##
## [[2]]
## [1] 1536 96
##
## [[3]]
## [1] 1536 48
##
## [[4]]
## [1] 1536 24
cov.error.ok.mesh <- matrix(NA, nrow = length(h_vector), ncol = length(beta_vector))
cov.error.folded <- matrix(NA, nrow = length(h_vector), ncol = length(beta_vector))
par(family = "Palatino")
layout(matrix(1:(length(beta_vector)*length(h_vector)),
nrow = length(beta_vector),
byrow = TRUE))
m_values <- c()
for (j in 1:length(beta_vector)) {
beta <- beta_vector[j]
alpha <- 2*beta
fract2beta <- alpha - floor(alpha)
nu <- alpha - 0.5
kappa <- sqrt(8*nu)/rho
tau <- sqrt(gamma(nu) / (sigma^2 * kappa^(2*nu) * (4*pi)^(1/2) * gamma(nu + 1/2))) #sigma = 1, d = 1
Sigma.folded <- gets_true_cov_mat(graph = graph.ok,
kappa = kappa,
tau = tau,
alpha = alpha,
n.overkill = n.overkill)
for (i in 1:length(h_vector)) {
h <- h_vector[i]
m <- min(20, ceiling((min(4*beta - 1/2,2) + 1/2)^2*log(h)^2/(4*pi^2*fract2beta)))
m_values <- c(m_values, m)
Sigma.ok.mesh <- matern.operators(alpha = alpha, # using overkill mesh
kappa = kappa,
tau = tau,
m = m,
graph = graph.ok,
type = type,
type_rational_approximation = type_rational_approximation)$covariance_mesh()
Sigma <- matern.operators(alpha = alpha, # this is the approximation
kappa = kappa,
tau = tau,
m = m,
graph = graphs[[i]],
type = type,
type_rational_approximation = type_rational_approximation)$covariance_mesh()
Sigma.approx <- A[[i]]%*%Sigma%*%t(A[[i]])
cov.error.ok.mesh[i,j] <- sqrt(as.double(t(graph.ok$mesh$weights)%*%(Sigma.ok.mesh - Sigma.approx)^2%*%graph.ok$mesh$weights))
cov.error.folded[i,j] <- sqrt(as.double(t(graph.ok$mesh$weights)%*%(Sigma.folded - Sigma.approx)^2%*%graph.ok$mesh$weights))
aux <- dim(Sigma.approx)[1]
reordering <- c(1,3:aux,2)
auxhalf <- reordering[length(reordering) %/% 2 + 1]
plot(Sigma.folded[auxhalf, reordering], col = "black", type = "l", main = bquote(beta == .(beta_label[j]) ~ "," ~ h == .(h_label[i])), xlab = "", ylab = "", lwd = 2, lty = 1)
lines(Sigma.ok.mesh[auxhalf, reordering], col = "red", lty = 2, lwd = 2)
lines(Sigma.approx[auxhalf, reordering], col = "blue", lty = 3, lwd = 2)
legend("topleft", lwd = 2, col = c("black", "red", "blue"), legend = c("true", "ok.mesh", "approx"), lty = c(1,2,3))
}
}## [1] 20 20 20 20 18 13 8 5 11 8 5 3 8 6 4 2 6 4 3 2 5 4 2 2 4
## [26] 3 2 1
## [1] 1.50 1.75 2.00 2.00 2.00 2.00 2.00
# folded
for (u in 1:length(beta_vector)) {print(coef(lm(log(cov.error.folded[,u]) ~ log(h_vector))))}## (Intercept) log(h_vector)
## 0.1731924 1.4498299
## (Intercept) log(h_vector)
## 0.4465996 1.6381700
## (Intercept) log(h_vector)
## 0.6274769 1.7777247
## (Intercept) log(h_vector)
## 0.6852249 1.8576150
## (Intercept) log(h_vector)
## 0.6937527 1.9022237
## (Intercept) log(h_vector)
## 0.6694985 1.9226518
## (Intercept) log(h_vector)
## 0.5039852 1.8960870
# ok.mesh
for (u in 1:length(beta_vector)) {print(coef(lm(log(cov.error.ok.mesh[,u]) ~ log(h_vector))))}## (Intercept) log(h_vector)
## 0.1592589 1.4441568
## (Intercept) log(h_vector)
## 0.4506415 1.6398464
## (Intercept) log(h_vector)
## 0.6679511 1.7933716
## (Intercept) log(h_vector)
## 0.7861186 1.8949556
## (Intercept) log(h_vector)
## 0.8361226 1.9559881
## (Intercept) log(h_vector)
## 0.8411164 1.9887793
## (Intercept) log(h_vector)
## 0.798770 1.998613
Press the Show button below to reveal the code.
loglog_line_equation <- function(x1, y1, slope) {
b <- log10(y1 / (x1 ^ slope))
function(x) {
(x ^ slope) * (10 ^ b)
}
}
guiding_lines <- matrix(NA, nrow = length(h_vector), ncol = length(beta_vector))
for (j in 1:length(beta_vector)) {
guiding_lines_aux <- matrix(NA, nrow = length(h_vector), ncol = length(h_vector))
for(k in 1:length(h_vector)){
point_x1 <- h_vector[k]
point_y1 <- cov.error.ok.mesh[k, j]
slope <- theoretical_rate[j]
line <- loglog_line_equation(x1 = point_x1, y1 = point_y1, slope = slope)
guiding_lines_aux[,k] <- line(h_vector)
}
guiding_lines[,j] <- apply(guiding_lines_aux, 1, mean)
}
# Generate default ggplot2 colors
default_colors <- scales::hue_pal()(ncol(guiding_lines))
# Create the plot_lines list with different colors for each line
plot_lines <- lapply(1:ncol(guiding_lines), function(i) {
geom_line(data = data.frame(x = h_vector, y = guiding_lines[, i]),
aes(x = x, y = y), color = default_colors[i], linetype = "dashed", show.legend = FALSE)
})
df <- as.data.frame(cbind(h_vector, cov.error.ok.mesh))
colnames(df) <- c("h_vector", beta_label)
df_melted <- melt(df, id.vars = "h_vector", variable.name = "column", value.name = "value")
p <- ggplot() +
geom_line(data = df_melted, aes(x = h_vector, y = value, color = column)) +
geom_point(data = df_melted, aes(x = h_vector, y = value, color = column)) +
plot_lines +
labs(title = "Covariance error using overkill mesh",
x = TeX("$\\hat{h}$"),
y = "Covariance Error",
color = expression(beta)) +
scale_x_log10(breaks = h_vector, labels = TeX(h_label_latex)) +
scale_y_log10() +
theme_minimal() +
theme(text = element_text(family = "Palatino"))Press the Show button below to reveal the code.
guiding_lines <- matrix(NA, nrow = length(h_vector), ncol = length(beta_vector))
for (j in 1:length(beta_vector)) {
guiding_lines_aux <- matrix(NA, nrow = length(h_vector), ncol = length(h_vector))
for(k in 1:length(h_vector)){
point_x1 <- h_vector[k]
point_y1 <- cov.error.folded[k, j]
slope <- theoretical_rate[j]
line <- loglog_line_equation(x1 = point_x1, y1 = point_y1, slope = slope)
guiding_lines_aux[,k] <- line(h_vector)
}
guiding_lines[,j] <- apply(guiding_lines_aux, 1, mean)
}
# Generate default ggplot2 colors
default_colors <- scales::hue_pal()(ncol(guiding_lines))
# Create the plot_lines list with different colors for each line
plot_lines <- lapply(1:ncol(guiding_lines), function(i) {
geom_line(data = data.frame(x = h_vector, y = guiding_lines[, i]),
aes(x = x, y = y), color = default_colors[i], linetype = "dashed", show.legend = FALSE)
})
df <- as.data.frame(cbind(h_vector, cov.error.folded))
colnames(df) <- c("h_vector", beta_label)
df_melted <- melt(df, id.vars = "h_vector", variable.name = "column", value.name = "value")
q <- ggplot() +
geom_line(data = df_melted, aes(x = h_vector, y = value, color = column)) +
geom_point(data = df_melted, aes(x = h_vector, y = value, color = column)) +
plot_lines +
labs(title = "Covariance error using folded",
x = TeX("$\\hat{h}$"),
y = "Covariance Error",
color = expression(beta)) +
scale_x_log10(breaks = h_vector, labels = TeX(h_label_latex)) +
scale_y_log10() +
theme_minimal() +
theme(text = element_text(family = "Palatino"))We used R version 4.4.0 (R Core Team 2024) and the following R packages: here v. 1.0.1 (Müller 2020), htmltools v. 0.5.8.1 (Cheng et al. 2024), INLA v. 24.6.27 (Rue, Martino, and Chopin 2009; Lindgren, Rue, and Lindström 2011; Martins et al. 2013; Lindgren and Rue 2015; De Coninck et al. 2016; Rue et al. 2017; Verbosio et al. 2017; Bakka et al. 2018; Kourounis, Fuchs, and Schenk 2018), inlabru v. 2.10.1.9010 (Yuan et al. 2017; Bachl et al. 2019), knitr v. 1.47 (Xie 2014, 2015, 2024), latex2exp v. 0.9.6 (Meschiari 2022), Matrix v. 1.6.5 (Bates, Maechler, and Jagan 2024), MetricGraph v. 1.3.0.9000 (Bolin, Simas, and Wallin 2023b, 2023a, 2023c, 2024; Bolin et al. 2023), patchwork v. 1.2.0 (Pedersen 2024), plotly v. 4.10.4 (Sievert 2020), plotrix v. 3.8.4 (J 2006), reshape2 v. 1.4.4 (Wickham 2007), rmarkdown v. 2.27 (Xie, Allaire, and Grolemund 2018; Xie, Dervieux, and Riederer 2020; Allaire et al. 2024), rSPDE v. 2.3.3.9000 (Bolin and Kirchner 2020; Bolin and Simas 2023; Bolin, Simas, and Xiong 2023), scales v. 1.3.0 (Wickham, Pedersen, and Seidel 2023), sf v. 1.0.16 (Pebesma 2018; Pebesma and Bivand 2023), tidyverse v. 2.0.0 (Wickham et al. 2019), xaringanExtra v. 0.8.0 (Aden-Buie and Warkentin 2024).